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A METHOD  FOR  THE  NUMERICAL  SOLUTION 
OF  A HEAT  CONDUCTION  PRCDLEM 


Prepared  by 
R.  2.  Eddy 


ABSTRACT:  A scheme,  due  to  J.  von  Nermann,  for  the  numerical  integration  of 

parabolic  or  hyperbolic  partial  differential  equations  is  here  set  up  for  use 
in  solving  a particular  problem  in  heat  conduction.  A crude  stability  analysis 
is  Given,  showing  that  the  effects  of  round-off  errors  tend  to  be  damped  out 
regardless  of  the  sise  of  the  time  step  — a situation  which  does  hot  hold 
true  for  some  better  known  integration  schemes. 
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This  report  contains  a description  of  a scheme  for  numerical  integration  of 
partial  differential,  equations  of  parabolic  or  of  hyperbolic  type,  due  to 
J.  von  Neumann.  Ibis  report  has  been  organized  under  Project  NR -044 -00 3, 
work  in  Numerical  Analysis,  sponsored  by  the  Office  of  Rival  Research. 
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A METHOD  FOR  THE  NUK"RICAL  SOLUTION 
OF  A HEAT  CONDUCTION  PROBLEM 


1 . Introduction 

In  the  course  of  development  of  a certain  high  velocity  projectile  it  became 
desirable  to  estia&te  the  temperature  of  the  projectile's  outer  shell  or  "skin/' 
during  the  course  of  flight.  One  way  of  accomplishing  this  was  to  obtain  a 
numerical  solution  of  the  partial  differential  equation  which  governs  heat  flow 
in  an  Infinite  plate , using  appropriate  boundary  conditions.  This  equation  vas 
non  linear  due  to  the  fact  that  the  thermal  properties  of  steel  vary  appreciably 
over  the  temperature  range  encountered.  At  the  inner  surface  the  boundary  con- 
dition was  taken  to  be  no  transfer  to  heat  across  the  surface,  while  at  the  outer 
surface  the  boundary  condition  took  into  account  heat  transfer  both  by  conduction 
to  or  from  the  boundary  layer  and  by  radiation  -nto  space.  (The  latter  effect 
turned  out  to  be  * lmost  negligible. ) 

A scheme  for  numerical  integration  suggested  by  Dr.  John  von  Neumann  of  the 
Institute  for  Advanced  Study  was  used  for  part  of  the  calculations  and  will  he 
described  here.  It  is  an  "implicit  method"  where  in  the  values  of  the  dependent 
variable  on  the  new  time  step  are  expressed  in  terms  of  each  other  and  must  be 
obtained  by  solving  a simple  system  of  linear  algebraic  equations;  its  chief 
interest  lies  in  its  stability:  round-off  errors  are  damped  out  regardless  of 

the  size  of  the  time  step. 

Numerical  results  which  have  been  obtained  in  connection  with  this  study 
can  be  provided  upon  request. 

2.  Equations 


The  governing  partial  differential  equation  was  taken  to  be 


ft  C ( T) 


with  ini tlul  condition 
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where  the  meaning  of  the  symbols  is  as  follows: 

T temperature  in  degrees  Rnnkine  (absolute  Fahrenheit) 
t time  in  hours 
x distance  in  feet 

K(T)  thermal  conductivity  in  BTU/hr.  ft.  °7 
C(T)  heat  capacity  in  BTU/lb.  °F 

f density  in  lb/ft3.  (assumed  constant) 
h(t)  coefficient  of  boat  transfer  between  gas  boundary  layer  and  steel  skin 
Git)  temperature  of  gas  boundary  layer 

Ta  effective  black  body  temperature  of  space,  assumed  for  simplicity  to  be 
constant  at  0°F  « 4j60°R. 
b net  coefficient  of  radiation  emlssivity 

The  function*  h(t)  and  G(t)  hod  been  calculated  previously  using  the  best 
data  available  for  a number  of  different  calculated  trajectories. 

The  functions  K(T)  and  C(T)  were  approximated  by  piece-wise  linear  functions 


like  those  in  Figure  1, 

K(T) 

fCtT) 

sv- 



fij  l'**  T 

Fij  L-b  T 

3.  Finite  difference  formulation. 

The  finite  difference  mesh  uaed  is  as  shown  in  Figure  2.  For  convenience  in 
satisfying  boundary  conditions,  the  mesh  is  laid  out  so  that  both  boundaries  fall 
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half  way  between  znesh  lines . The  number  of  vertical  mesh  lines  (size  of  Ax) 
was  determined  largely  by  the  memory  capacity  of  the  oc  touting  machine  used! 
the  Mark  II  Relay  Calculator  at  the  Naval  Proving  Ground,  Dahlgren.  Virginia. 
A superscript  on  a symbol  refers  to  tho  time  coordinate,  while  a subscript 
refers  to  the  space  coordinate.  Thus 

r‘  = T(0-Y>*&x,  tat) 

•j-'fX.  v 


The  presence  of  non-integer  indices  is  due  partly  to  the  fact  th«?t  the  first  mesh 
line  in  from  the  left  boundary  (Figure  2)  lies  at  x “ 1/2  Ax  instead  of  at  x =Ax, 
and  partly  to  the  extensive  use  of  averaging  or  forming  linear  combinations  of  the 
values  of  some  variable  at  several  mesh  points  in  order  to  define  the  value  of 
another  variable  at  an  Intermediate  point.  The  reason  for  the  letter  la,  of  course, 
to  Secure  greater  accuracy. 

The  following  auxiliary  quantities  are  introduced! 
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With 

these  definitions  one  has 
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r ccr)k%‘  = rc(t?j  { t-:::  - }M 

Substitution  of  those  into  (l)  leads  to  the  following  system  of  equations 

c t"a.  , „ t*Wt 


(15) 


where 


'V'-Va  /V-Vk  - 0/-4fe  /»</*  ~ 4/-J/a 

a.*>  = (%-C + uj  ) 

6j-x,  z Uj//y£  ( (J*  2,3,*,s) 

n.j,x  -•  WZ  J 


(16) 


(17) 


Prom  (3)  and  (13)  one  has  7>«/t  - %a  ®o  that  (16)  allows  ''Jf/i  to  be 

expressed  in  the  for*  J7/^  r ^MVx  yVt  - 1°  general 


‘ i 7l/z  ~ V- 
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where 


| (J '«>■*, 


3>  (19) 


idien 
obtained 
gives 


•fg  ~ /^VVv, 

Vj-Xx.-  AjSIx  yj.’/x  ~ Vjt'tx  I’fj-J/x.  j 

J » 1 this  recursion  schema  breaks  down  sinoe  ▼_  ■ 0)  however  the  equation 
ned  from  (16)  by  nultiplying  through  by  ▼j„i>  then  putting  j = 1 and  v“  « 0, 


(20) 


Elimination  of  Jj/i  and  7'A  by  means  of  (18)  with  J ■ 3 and  ] ■ 2,  respec- 
tively, leaves  a single  equation  in  ^Vt  *61  ch  can  be  solved  inmediately.  Sub- 
stitution of  this  value  into  the  expressions  already  set  up  from  (18)  then  gives 
the  values  of  the  remaining  7/-**  for  J * 5,  4,  3,  2.  Prom  this  the  desired 

temperatures  are  then  found  with  the  aid  of  (6) t 


'j-*A  °*  /'-- 


v-vi 


(21) 


4.  Stability  Argument 

As  mentioned  in  the  introduction,  one  of  the  interesting  features  of  this 
integration  scheme  is  that  errors  due  to  round-off  are  damped  out  no  matter  how 
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large  /It  ie  taken  — soaethJ  g which  lc  not  true  of  some  simpler  schemes.  The 
following  crude  argument,  although  not  a rigorous  proof,  is  of  a kind  which  has 
proved  to  be  quite  useful  and  fairly  accurate  in  practice. 

The  first  step  is  to  ignore  the  variability  of  the  coefficients  of  (l);  this 
leads  to  a partial  difference  equation  with  constant  coefficients  and  homogeneous 
boundary  conditions,  which  is  satisfied  by  any  round-off  error  function,  and 
which  can  be  solved  by  separation  of  variables  t 
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■'ft. 


- O 


(22) 


where 

r _ K 

-Jf  C-fA*;* 

(23) 

Let 

tA=  r 

(24) 

where  U is  a funotion  of  x alone  and  V is  a function  of  t alone. 
(24)  into  (22)  and  writing  (A-OAA+O  for  the  separation  constant, 

Substituting 

one  has 

(25) 

-2(l-Jr)0j  + Li‘.,'0  . 

(26) 

To  solve  (26) , put 

1 

-A  = f 

(27) 

Uj  = ^ 

(28) 

which  leads  to 

^ - eT 

(29) 

so  that 

\Jy  r <*  COJ  Jf  + fl  Si*  Jf  . 

(30) 

The  boundary  conditions  then  determine  a set  of  values  of  such  that  (30)  is 

an  acceptable  solution  of  (26).  These,  in  turn,  determine  the  eigenvalues 


(31) 


which  are  never  negative,  since  r is  positive. 


i A < i 


Therefore 


(32) 
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for  ell  values  of  A encountered  so  that  the  solution  of  (25) 

VC*V0(-H7-y  (33, 

regains  bounded  as  t increases.  This  means  that  the  finite  difference  scheme 
considered  here  i6  always  stable,  regardless  of  the  else  of  ^t. 
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